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abstract 

We present a global magnetohydrodynamic (MHD) three dimensional simulation of a nonradiative accretion 
flow originating in a pressure supported torus. The evolution is controlled by the magnetorotational instability 
which produces turbulence. The flow forms a nearly Keplerian disk. The total pressure scale height in this disk 
is comparable to the vertical size of the initial torus. Gas pressure dominates only near the equator; magnetic 
pressure is more important in the surrounding atmosphere. A magnetically dominated bound outflow is driven 
from the disk. The accretion rate through the disk exceeds the final rate into the hole, and a hot torus forms inside 
10 rg. Hot gas, pushed up against the centrifugal barrier and confined by magnetic pressure, is ejected in a narrow, 
unbound, conical outflow. The dynamics are controlled by magnetic turbulence, not thermal convection, and a 
hydrodynamic a model is inadequate to describe the flow. The limitations of two dimensional MHD simulations 
are also discussed. 

Subject headings: accretion — accretion disks — instabilities — MHD — black hole physics — X-ray: stars - 
binaries: 
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1. INTRODUCTION 

Observations of underluminous accretion-powered X-ray 
sources pose a difficulty for accretion theory. The Galactic cen- 
ter is a particularly compelling example: despite the inferred 
presence of a 2 million solar mass black hole, and the avail- 
ability of a substantial gas reservoir, only a small fraction of 
the expected radiation is detected. In a standard Keplerian ac- 
cretion disk, liberated mechanical energy is promptly radiated, 
and there is a direct relationship between the accretion rate M 
and the luminosity L. A number of alternatives to the standard 
disk have been put forward, which are classified by the ultimate 
repository of the orbital energy. In Advection Dominated Ac- 
cretion Flows (ADAFs, e.g., Narayan & Yi 1994), the energy is 
retained by the plasma and advected into the hole. In contrast, 
Blandford & Begelman (1999) proposed that the bulk of the ac- 
creting gas and energy might be carried off by a wind. A more 
recent variation, referred to as Convective Dominated Accre- 
tion Flow (CDAF; e.g., Quataert & Gruzinov 2000), brings en- 
hanced viscous heating to the fore, leading to a strong convec- 
tive flow that stifles the accretion. Since the common feature of 
all these models is negligible radiative losses, we shall refer to 
them generically as Non-Radiative Accretion Flows (NRAFs). 

Earlier simulations of NRAFs have used several simplify- 
ing approximations. Most (e.g., Igumenshchev & Abramowicz 
1999, 2000, hereafter IA99, lAOO; Stone, Pringle, & Begel- 
man 1999, hereafter SPB) have assumed the flow evolves due 
to a fortified kinematic viscosity, v. While illuminating some 
aspects of NRAFs, there are substantial drawbacks to this ap- 
proach. Most importantly, the results strongly depend upon the 
specific recipe adopted for v. The origin of angular momen- 
tum transport is not unknown, however. It arises from the mag- 
netorotational instability (MRI; Balbus & Hawley 1991), and 
hydrodynamic and magnetohydrodynamic (MHD) flows have 
fundamentally different properties. 

Stone & Pringle (2001; hereafter SP) include MHD in their 
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axisymmetric simulations. The stress required to drive the ac- 
cretion flow emerges self-consistently; it is not a free parameter 
But the two-dimensional (2D) restriction is a significant limita- 
tion. First, the anti-dynamo theorem (e.g., Moffat 1978) pre- 
vents the indefinite maintenance of the poloidal magnetic field 
in the face of dissipation. Second, axisymmetric simulations 
tend to over-emphasize the "channel" mode (Hawley & Bal- 
bus 1992) which produces coherent internal magnetized flows 
rather than the more generic MHD turbulence. Consequently, 
a fully self-consistent NRAF model requires three-dimensional 
(3D) MHD. In this Letter we present and discuss the results of 
a prototype 3D NRAF simulation. 

2. THE SIMULATION 

The simulation evolves the equations of ideal MHD, i.e., the 
continuity equation, the induction equation, the Euler equation, 

p| + (Pv.V)v=-v(p+c1 + ^)-PV<I>+(A.v)b, 

(1) 



and an internal energy equation 
3pe 



-^ + V.(pev) = -(P+Q,)V-v. 
of 



(2) 



The variables have their usual meanings; Ci is an explicit artifi- 
cial viscosity (Stone & Norman 1992a), and <l> — —GMj (r — rg) 
is the Paczynski & Wiita (1980) pseudo-Newtonian potential. 
We use units with GM = rg = I. The equation of state is adi- 
abatic, P = pe(r — 1), with F = 5/3. Radiation transport and 
losses are omitted. Since there is no explicit resistivity or shear 
viscosity, the gas can heat only by adiabatic compression, or by 
the artificial viscosity Q,. We evolve the equations using time- 
explicit Eulerian finite differencing with the ZEUS algorithms 
(Stone & Norman 1992a; 1992b; Hawley & Stone 1995; Haw- 
ley 2000). 
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The major difficuhies of a 3D simulation include the large 
number of grid zones needed to resolve an extended spatial do- 
main, and the number of time steps required to evolve the prob- 
lem over several orbital periods at large radius. In this simula- 
tion, we use a somewhat restricted resolution of 128 x 32 x 128 
grid zones in cylindrical coordinates (/?,(]), Z). The radial grid 
extends from R — 1.5 to R = 170. There are 36 equally-spaced 
zones inside R = 15, and 92 zones outside this point that in- 
crease logarithmically in size. Similarly, 50 Z zones are equally 
spaced between —10 and 10 with the remainder of the zones 
logarithmically stretched to the Z boundaries at ±60. The az- 
imuthal domain is limited to 7t/2. Hawley (2001) found that 
reduced angular coverage preserves qualitative features of a 
simulation with only a small reduction in magnetic energy and 
stress. The radial and vertical boundary conditions are sim- 
ple zero-gradient outflow conditions; no flow into the compu- 
tational domain is permitted. The (|) boundary is periodic. The 
magnetic field boundary condition is set by requiring the trans- 
verse components of the field to be zero outside the computa- 
tional domain, while the perpendicular component satisfies the 
divergence-free constraint. 

The prototype simulation starts with of a constant specific 
angular momentum (£) torus with a pressure maximum at 
r = 100 (similar to model F of SP) and an inner edge at 
r = 75. The initial magnetic field consists of poloidal loops 
lying along isodensity contours with a volume-averaged en- 
ergy of p = Pgas/Pmag = 200. The average specific energy of 
the magnetofluid in the torus is —4.6 x 10^^ (in units where 
GM = = 1), i.e., it is bound. The simulation is run out to 5 or- 
bits at the initial pressure maximum. This corresponds to 1429 
orbits at the marginally stable orbit, 7? = 3. The simulation re- 
quires over 1.7 million timesteps with an average Af = 0.0152. 

3. RESULTS 

During the first orbit the field in the torus is amplified by the 
MRI and by shear, causing the torus to expand. The MRI acts 
most visibly near the equatorial plane, where the field is pre- 
dominantly vertical, and the long-wavelength, nearly axisym- 
metric modes of the MRI grow rapidly. The linear growth phase 
ends shortly after one orbit, and the magnetic energy has in- 
creased to p « 2-10. 

Over the next orbit, low-m spiral arms of gas accrete from 
the inner edge of the torus, forming a vertically thin and very 
nonuniform disk. Strong magnetic fields surround the gas. Due 
to the initial field topology, a current sheet forms near the equa- 
tor which proves unstable to vertical oscillations. As more gas 
leaves the torus, the disk fills out and thickens. Low density 
material is driven off the forming disk to create a backflow. 

As time advances, the initial constant-^ torus evolves to a 
nearly Keplerian angular momentum distribution. There is net 
accretion inside ofR= 100, and net outflow beyond this point. 
Inside of R = 100 the NRAF forms a modestly thick, nearly 
Keplerian disk. In the disk, accretion results from MHD turbu- 
lence, but the final accretion rate into the black hole need not 
match the supply rate precisely, and it does not. Because the 
inflow is hot, geometrically thick, and nearly Keplerian, inflow 
is difficult past the marginally stable orbit. To accrete the gas 
must pass through a narrow gap in the centrifugal barrier at the 
equator. Gas splashes off the centrifugal barrier near the hole, 
forming a pressure-supported torus and creating a backflow that 
adds to a growing low density envelope around a higher density 
Keplerian core. 



Between orbits 4 and 5, the disk is well-formed and reason- 
ably steady. Figure 1 is a series of {R,z) contour plots of az- 
imuthally averaged values &Xt = 5 orbits showing that most of 
the gas is located in a modestly thick disk surrounded by a low 
density, highly magnetized atmosphere. Near the equator gas 
pressure is dominant, i.e., p > 1, but in the surrounding region 
P < 1. The gas pressure scale height is // « 7-10, decreasing 
rapidly inside R = 10. The total pressure (gas plus magnetic) 
is much smoother, and has a scale height H « 20 at /? = 100 
that decreases slowly inward. The gas density, vertically aver- 
aged over one gas pressure scale height, is nearly constant from 
R = 50 down to R^ 10, where it increases rapidly to a peak at 
R = 5.5. Inside /?= 10, the disk resembles a pressure supported 
torus. 

The specific angular momentum 1{R) tracks a Keplerian dis- 
tribution throughout the disk, although in places it is up to 10% 
below the Keplerian value. The angular velocity Q. is constant 
along cyhnders through the main portion of the disk, consistent 
with P = P{p) there. In hydrodynamic simulations, surfaces of 
constant entropy S coincide with those of specific angular mo- 
mentum t. This corresponds to marginal stabihty by one of the 
H0iland criteria. As already noted by SP, with MHD the level 
surfaces of S and t do not line up, particularly inside the main 
disk. 

Magnetic fields in the disk produce Maxwell stress, MHD 
turbulence, and angular momentum transport. Between orbit 4 
and 5 the ratio of the vertical- and (|)-averaged Maxwell stress 
to the gas pressure ranges from 0.1 to 0.2 inside of R = 100. 

The mass inflow and outflow rates through every cylindrical 
radius are computed as a function of time. An average over 
the last orbit shows that these two are nearly equal, although 
inflow exceeds outflow (Fig. 2). Of course the instantaneous 
inflow and outflow rates merely represent the flow produced 
by the MHD turbulent velocity fluctuations. The net accretion 
rate, M, is due to the resulting drift velocity, which is always 
considerably smaller. 

At the end of the run the total mass on the grid has decreased 
by 5.4%. Roughly 51% of this leaves through the outer radial 
boundary, 43% through the upper and lower Z boundaries, and 
the remaining 6% is accreted into the black hole. A higher res- 
olution torus simulation (Hawley & Krolik 2001) shows that 
the magnetic stress can remain large down to and beyond the 
marginally stable orbit. This effect could increase M, but the 
present simulation is not sufficiently well-resolved to address 
this point. The flow through the outer radial boundary is a con- 
sequence of the gain in t in the outer part of the torus. The 
flow through the Z boundaries is driven from the accretion disk 
by gas and magnetic pressure. Part of this is an unbound, high 
temperature hollow conical outflow confined to the axis region 
by surrounding magnetic pressure. This can be seen in the gas 
pressure plot of Figure 1 . The remainder of the outflow remains 
bound, although its specific energy is greater than that of the 
initial torus. 

The initial average specific energy for the gas is —4.6 x 10~^. 
At 5 orbits this has decreased to —4.9 x 10~^ as the gas re- 
maining on the grid becomes more bound. The total thermal, 
magnetic, poloidal kinetic, and orbital energies of the gas that 
remains on the grid all increase with time. The magnetic energy 
has increased the most, receiving 42% of the total net energy in- 
crease; 95% of this is in toroidal field. The thermal, orbital and 
kinetic energies comprise 25%, 21%, and 12% of the increase 
respectively. 
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4. DISCUSSION 
4.1. Comparison with a models 

Full 3D MHD simulations are difficult and expensive; are 
they necessary? With respect to the hydrodynamical alternative 
the answer is clearly yes. Such simulations (IA99; SPB; lAOO) 
have shown that the results depend strongly on the magnitude 
of a and the form assumed for the stress. Thus, a self-consistent 
stress is essential for understanding NRAF dynamics. 

One should not lose sight of the fact that this is a high 
Reynolds number turbulent system, not a low Reynolds num- 
ber laminar flow. For example, the result that large Shakura- 
Sunyaev a (v = occ^ /Q.) accretion flows are stable, laminar and 
accrete fully into the central hole (IA99; lAOO) is a Uteral con- 
sequence of using a Navier-Stokes viscosity for the stress. With 
stress from MHD turbulence a ~ 1 necessarily imphes velocity 
fluctuations with speeds comparable to the sound speed Cs on 
scales of order the size of the system. 

A second hydrodynamic result is that NRAFs are driven to 
a state of marginally convective instability (SPB; Quataert & 
Gruzinov 2000; lAOO). This is the basis for CDAF solutions. 
But this cannot hold in MHD. Recall the form of the H0iland 
criteria that appUes to adiabatic perturbations in magnetized 
flow (Balbus 1995): 
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>0 
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dR 



(3) 



>0 (4) 



where ^ -(3/5p)VP- Vlnfp^VS is the Brunt- Vaisala fre- 
quency. Note the symmetry between the angular velocity gradi- 
ents and the entropy gradients, and the absence of any stabiliz- 
ing epicyclic frequency. At least one of these criteria must al- 
ways be strongly, not marginally, violated in an accretion flow: 
this is what is responsible for the existence of the turbulent 
stress. Convection does not arise in addition to this fundamental 
instability, it is part of the instability that produces the anoma- 
lous stress in the first place. In particular, one cannot argue for 
the marginal stability of a rotationally modified convective pro- 
cess independently of what is called a. Whatever the sign of the 
Brunt-Vaisala frequency, the instability is essentially the same 
one at work in a standard Keplerian disk. Indeed, in contrast to 
hydrodynamical treatments, we find that the MHD flows show 
no tendency at all toward marginal stability. 

The 3D MHD simulation does show large scale turbulent 
rolls and outflows quahtatively similar to those reported in 
small a hydrodynamic simulations. In the hydrodynamic simu- 
lations, however, all the rotational energy released by the stress 
goes immediately into heat. Some of this heat returns to kinetic 
energy in the form of large-scale convection. This implies that 
large-scale shear is converted into heat, producing large-scale 
convective rolls, which then prevent the accretion process that 
was driving the heating in the first place. In MHD, the orbital 
energy goes directly into turbulent kinetic and magnetic ener- 
gies. The magnetic energy constitutes a particularly important 
component of the total energy budget. Heating results from dis- 
sipation at small scales. In a real plasma, this will depend on 
microscopic resistivity and viscosity. In the simulation, turbu- 
lent rolls are not thermal in origin. More generally the kinetic 
energy of convective motions is small compared with the dy- 
namical activity associated with MHD instabilities and with the 
large scale flow itself. 



We also remark that in any turbulent flow, the net accretion 
rate will always be small compared to any locally sampled in- 
stantaneous value. As Figure 2 shows, inward and outward ac- 
cretion rates almost exactly cancel in our purely dynamical sim- 
ulation. This is not a unique property of inwardly-transported 
convective energy fluxes. It is a hallmark of turbulence. In both 
viscous hydrodynamical and MHD treatments, turbulent eddy 
velocities are much larger than the residual inward drift veloc- 
ity. It is precisely this ordering that allows the a formalism for 
the mean flow dynamics to be derived from the equations of 
motion themselves (Balbus & Papaloizou 1999). 



4.2. Comparison with 2D MHD models 

Next we compare the present simulation with the 2D models 
of SP. In the equivalent run (SP's Run F), the initial infaU phase 
during orbits 2-3 was relatively smooth and dominated by the 
channel flow of the MRI. Subsequent evolution was turbulent. 
In 3D, the initial smooth phase does not occur; the evolution is 
always turbulent. The 3D flow also generates greater outflow 
from the disk. 

After 3 orbits in both 2D and 3D, the inner regions of the flow 
are dominated by MHD turbulence driven by the MRI. Inter- 
estingly, there appear to be few differences between time- and 
angle-averaged quantities in either case. For example, com- 
parison of the contours in Fig. 1 with Figure 4 in SP show 
that both produce a nearly barytropic disk near the midplane 
with contours parallel to cylindrical radii, but with no cor- 
relation between the specific angular momentum and entropy. 
Moreover, the net M values are similar: 3 x 10~^ in 2D versus 
^ 2 X 10^^ in 3D, in units of the initial mass of the torus per 
orbit at the initial pressure maximum. Radial profiles of time- 
averaged data in the main disk show many similarities: in both 
cases Pgas ~ lOPmag and V(^>Cs> VAifvm at aU radii (see figure 
6 in SP). The most obvious differences between the 2D and 3D 
models are near the axis, and these are attributable to spherical 
2D versus cylindrical 3D inner boundary geometries. 

Eventually all axisymmetric models are limited by the an- 
tidynamo theorem. Toward the end of the SP simulations the 
turbulence in the initial torus dies down, and remains only in 
the flow close to the black hole. Higher resolution is possible 
in 2D, but this only prolongs the duration of the turbulent phase 
before eventual decline. In addition, the toroidal field MRI can- 
not be simulated in axisymmetry, which limits 2D models to 
initial conditions that contain large-scale poloidal field. 

5. CONCLUSION 

We have performed a three-dimensional MHD simulation of 
an NRAF originating from a constant-^ torus located at 100 
gravitational radii. The resulting flow does not resemble a clas- 
sical ADAF. Most of the infalhng mass and released energy 
does not end in the hole. In this sense the flow might be bet- 
ter described as an inflow/outflow solution of the type outlined 
by Blandford & Begelman (1999). The general appearance of 
the flow is similar to the low-a hydrodynamic models of SPB, 
IA99 and lAOO. Of course this is what one should expect for 
any hot, rotating, nearly Keplerian accretion flow. We have ar- 
gued that the CDAF picture is also inappropriate because the 
the classic H0illand criteria do not apply to an MHD flow. The 
outflow is not driven by thermal convection, but by the action 
of the same instability that accounts for a, namely the MRI. 
Marginally stabihty is not attained. 
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We conclude by considering some directions for future work. 
In the present simulation an internal energy equation is used and 
the only source of dissipative heating is the artificial viscosity 
Q,. Some amount of the total energy is inevitably lost in this 
formulation. A portion of this can be recaptured through arti- 
ficial resistivity, although this made no significant difference in 
flow structure or accretion rates when tested in 2D (SP). In any 
case, a more detailed study of the energy flow in a 3D simula- 
tion would be enlightening. Higher resolution is needed in the 
inner region to capture the details of the inner disk boundary 
and the flow past the marginally stable orbit. Models with a 
greater radial domain would allow the inflow to proceed over 
several decades in distance and to liberate greater amounts of 
gravitational energy. High-resolution axisynmietric MHD sim- 



ulations can play an important role here in mapping out the 
types of flows that are possible. Longer evolutions would pro- 
vide greater confidence that the results have lost memory of the 
initial conditions, but, because of the antidynamo theorem, true 
steady states can orfly be obtained in 3D. 
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Fig. 1 . — Azimuthally-averaged contour plots of indicated quantities in the iimer region of the computation at f = 5 orbits. Plots are gas density, log of the gas 
pressure, log of the magnetic pressure, log of the total pressure, log of the angular velocity D., specific angular momentum £, the entropy, S = \d(P/p^/^), and the 
Maxwell stress, Tr^ = —BrB^/Ati. 
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Fig. 2. — Amplitude of the total mass inflow rate (solid line), mass outflow rate (dotted line), and net accretion rate (dashed line) through cylindrical radius R, 
averaged in time between orbit 4 and 5. 



